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Abstract. - Cosmological N-body simulations aim to calculate the non-linear gravitational 
growth of structures via particle dynamics. A crucial problem concerns the setting-up of the 
initial particle distribution, as standard theories of galaxy formation predict the properties of 
the initial continuous density field with small amplitude correlated Gaussian fluctuations. The 
discretisation of such a field is a complex issue and particle fluctuations are especially relevant 
at small scales where non-linear dynamics firstly takes place. In general, most of the procedures 
which may discretise a continuous field, gives rise to Poisson noise, which would then dominate 
the non-linear small-scales dynamics due to nearest-neighbours interactions. In order to avoid 
such a noise, and to consider the dynamics as due only to large scale (smooth) fluctuations, 
an ad-hoc method (lattice or glassy system plus correlated displacements) has been introduced 
and used in cosmological simulations. We show that such a method gives rise to a particle 
distribution which does not have any of the correlation properties of the theoretical continuous 
density field. This is because discreteness effects, different from Poisson noise but nevertheless 
very important, determine particle fluctuations at any scale, making it completely different from 
the original continuous field. We conclude that discreteness effects play a central role in the 
non-linear evolution of N-body simulations. 



The purpose of cosmological N-body simulations is to calculate the non-linear growth of 
structures in the universe by following individual particles trajectories under the action of 
their mutual gravity |l], |^, B|. These particles are not galaxies but are meant to represent 
some sorts of collisionless clouds of elementary dark matter particles. In order to make them 
move, one must calculate the force acting on each of them due to all the others. In general 
one may find several algorithms which speed up the N 2 sum necessary to compute the force 
on each particle. In this paper we study simulations from the Virgo project || which are done 
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with an "adaptative P 3 M" scheme This scheme is a combination of a PM (i.e. meshed 

based scheme) and a PP (i.e. particle-particle based scheme). Briefly, this consists of assigning 
the particle mass on a mesh so as to obtain a density held on the mesh itself. One can then 
calculate the gravitational potential on that mesh by solving the Poisson equation, and finally 
the potential at the location of a particle is determined by an interpolation of the values on the 
mesh. For a better resolution, a true calculation of the force due to nearby particles is done. 
In highly clustered regions, since this last calculation can be quite long, additional higher 
resolution meshes are placed. The force used is not a pure r~ 2 one: instead one smoothes 
it at small r by choosing for instance a force proportional to (r 2 + e 2 ) -1 in order to avoid 
"collisions" between close particles. This brings us to an important hypothesis sometimes 
made in cosmological N-body simulations: with a softened force and a proper choice of the 
softening parameter e, the evolution of the N-bodies should be the same as the evolution of a 
continuous density held (made of a huge number of particles behaving like a fluid) under its 
own gravity. With this in mind, one interprets collisions, or more correctly strong scattering, 
between nearby particles in the simulation as due to the discretisation of the density held and 
therefore artificial || However, it is important to note that, in order to satisfy the above 
hypothesis, e should be at least as large as the mean inter-particle distance (A) || this 
is not always the case as for instance in the Virgo project where e = 0.036Mpc/ft. [H] and 
(A) w Q.hMpc/h (see below) . 

The philosophy behind cosmological N-body simulations is to reproduce galaxy structures 
that we see today through redshift surveys by choosing some good parameters (cosmological 
and numerical) for the evolution of the system and especially fine-tuned initial conditions 
(properly normalised to the primordial fluctuations) given by some theories. This last point 
is important because one is not interested in some general asymptotic behaviors or quasi 
stationary states independent of the initial conditions, but in the state of the system after a 
relatively short time compared with the dynamical time of the simulation. 

Initial conditions (hereafter IC) are created from theories typically based on inflation and 
the subsequent evolution of matter, which are able to predict the properties of a continuous 
density field with correlated fluctuations || . Such a density field is usually Gaussian, and hence 
all the statistical properties are contained in the two-point correlation function (hereafter CF) 
or its Fourier transform, the power spectrum (hereafter PS). These two functions depend on 
the kind of dark matter which is the relevant one on large scales: hot (relativistic), cold 
(non-relativistic) or warm (a mixture of the two). However the main point for what concerns 
us is that a continuous and smooth density field with correlated density fluctuations is given 
as IC. If one wants to study the time evolution of this field with N-body simulations based 
on particle dynamics, it is then necessary to discretise the field. This means that one has 
to create a particle distribution which is representative of the continuous density field and to 
control any finite size effect. Note that a crucial point with respect to structures formation, 
is that non-linearities firstly develop at small scales where discreteness effects are important. 
Firstly, we briefly explain the standard way of carrying out the discretisation and then we 
analyse some of the real space statistical properties of these distributions, restricting ourself 
to a Virgo Standard Cold Dark Matter (hereafter CDM) simulation with 256 3 particles 
The questions we address are: In which range of scales the initial particle distribution used 
in the standard cosmological simulations is representative of CDM-like density fluctuations? 
What drives the non-linear dynamics of structures formation: the discreteness effects due to 
the short NN interaction or the large scale (continuous) distribution of density fluctuations? 

The standard ad- hoc procedure for setting up IC is described in |^, ||, H]. For galaxy 
structures formation problems, the IC generation splits into two parts. The first is to set 
up a "uniform" distribution of particles, which can represent the unperturbed universe. The 
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second is to impose density fluctuations with the desired characteristics. In this context, one 
faces the general problem of how to discretise a continuous density field. There are different 
procedures (e.g. random sampling, threshold sampling, etc.) which can be chosen and they 
result in different point distributions. Clearly one should have some physical reasons to choose 
one or another since any procedure introduces some discreteness effects, like Poisson noise, 
which could then play the main role in the non-linear dynamics of the system. For instance, 
in a Poisson distribution the dominant part of the gravitational force acting on an average 
particle is due to its NN Jl(| 111], [l2| . If a simulation is run from such IC the fluctuations grow 
rapidly into non-linear objects at small scales. This happens because a Poisson distribution 
is statistically isotropic only on scales larger than the average distance between NN jio] 11 



and the three-point correlation function is non-zero at these small scales: such a situation 
avoids the perfect cancelation of the net gravitational force due to close particles. Instead, 
in cosmological N-body simulations, one would like to simulate a system where the main 
contribution for non-linear structures formation problem, is due to the large scale distribution 
of the other particles and not to local NN interactions Q . 

To overcome the fact that a Poisson distribution is only "statistically" isotropic the most 
widely used solution to this problem has been to represent the unperturbed universe by a 
regular cubic grid of particles [fil 0, An infinite lattice, or a lattice with periodic boundary 
conditions, is "gravitationally stable" because of symmetry reasons. However a lattice is a 
distribution with fluctuations at all scales and non-trivial correlations. One can show that 
the unconditional variance in spheres decays as er 2 (r) = (AN(r) 2 ) / (N(r)) 2 ~ r~ 4 (where 
(N(r)) ~ r 3 is the average number of points in a randomly chosen sphere of radius r [p^[) and 
this is due to the fact that a lattice is a strongly ordered and correlated system at all scales. 
The two-point CF is such that £(fi,7^) = £{r\ — r 2 ) ^ £(|ri — r 2 |) because it is not invariant 
for space rotation |H: A lattice breaks space isotropy. Moreover, the grid-like system has the 
disadvantage to introduce a strong characteristic length on small scales - the grid spacing - 
and it leads to strongly preferred directions on all scales. 

For the latter reasons, a second way to generate an "uniform background" is by means of the 
following procedure. One starts from a Poisson distribution and then the N-body integrator is 
used with a repulsive gravitational force in such a way that, after the simulation is evolved for 
a sufficiently long time, the particles settle down to a glass-like configuration in which the force 
on each particle is very close to zero, i.e. a global stable configuration is found which has no 
preferred directions . The resulting distribution is very isotropic but it is still characterized 
by long-range order of the same kind of a lattice. As for the lattice the main characteristic 
of such a distribution is the presence of an excluded volume: two particles cannot lie at a 
distance smaller than a certain fixed length scale |jL4l. In the lattice this scale is the grid 
space, for a glass such a distance depends on the number of points one has distributed in a 
given volume. The fact that a lattice is ordered is due to the existence of the deterministic 
small scale distance. The unconditional variance scales as a 2 (r) ~ r~ a where 3 < a < 4, 
and it is again a strongly correlated system [Q. Its two-point CF depends on the detailed 
procedures used to generate the glass distribution. However it is possible to show that this 
class of distribution has P(k) ~ k a (with a > 0) for k — > and hence P(0) = like the 
Harrison- Zeldovich PS Q. 

Given a "suitably unperturbed" particle distribution, any desired linear fluctuation distri- 
bution can be in principle generated using Zeldovich approximation |Q . Basically, one assumes 
that the initial background is uniform with average density p(r) = po, without fluctuations, 
and then one applies a displacements field u(r). Because of the conservation of total mass 



4 



EUROPHYSICS LETTERS 



after the displacements, one can apply a continuity equation which gives 

p(r) - po + V • (potT(f)) = . (1) 

If statistical homogeneity and isotropy of the displacements field u(x) is assumed, then on 
going to Fourier space and considering the expectation value of the square modulus of the 
density fluctuations (the PS), EqJ| leads to P(k) = = k 2 (\u(k)\ 2 ) = k 2 P u {k). As any PS 
cannot diverge faster than fc~ 3 for k —> 0, we have that lim^^o P(k) ~ k n with n > —I. 
This means that, for instance, one cannot create a density fluctuation field with P(k) which 
diverges as k n and n < — 1, or with CF which goes as r~ a for r — > oo and a < 2 only by 
applying a displacement field to a uniform background. Moreover one does not start from 
a continuous field with no fluctuations, but with a particle distribution, which has its own 
fluctuations and correlations (which in the lattice's and glass's cases are long-range). In such a 
situation one has to ensure that the correlations among density fluctuations implemented by the 
displacements field are larger than the intrinsic fluctuations of the original particle distribution, 
and that the large-scale fluctuations dominate the non-linear small scales clustering instead of 
nearest-particle interactions. Only if one considers very large scales, and/or a displacements 
field which introduces correlations which are larger than the intrinsic one of the original 
distribution, one can recover the PS as in Eq.Q. Otherwise in a certain range of small enough 
scales, the point distribution is dominated by discreteness effects, which in this context can be 
seen as finite-size effects and which are important for what concerns the small-scale non-linear 
structures formation. 

In order to clarify this question, we have checked numerically the kind of correlations 
introduced in the system with the Zeldovich approximation: We have analysed a CDM 
simulation with Qq = 1 || as an example, but our result are generally valid for any other 
particular model chosen, as they involve the same method for setting-up initial conditions. 
Note that we have analysed the statistical quantities for the entire simulation (N = 256 3 
particles) and not for a sub-sample of it, in order to avoid introducing any kind of sampling 
noise and we have used periodic boundary conditions. The pre-initial particle distribution 
(which would represent the "uniform background" as previously discussed) consists of replicas 
of a N = I0 6 particles glass distribution generated by the procedure discussed above. As they 
are generated with periodic boundary conditions they can be tiled to make larger distributions. 
After the Zeldovich displacements one obtains the IC which correspond to a redshift z = 50. 
(In general the particles are very little displaced with respect to their initial positions: that is 
the average distance between NN remains unperturbed.) 

It is worth noticing that in general |IJ ||] the correlation properties of the initial particle 
distribution of N-body simulations have been discussed through the analysis of the PS of 
the density fluctuations. The PS is the Fourier transform of the real space CF £(r) and 
hence it contains the same information, clearly if finite size effects and statistical noise have 
been properly taken into account. For the comparison with theoretical CDM models, which 
usually give the PS of density fluctuations ||, we have computed the real space properties: 
f (r) = J °° P(k) s -i^ k 2 dk and a 2 {R) = ^ f™ W(k, R) 2 P{k) d 3 k where W(k, R) is the 
Fourier transform of the sample's window function (a real space sphere in our case - which is 
defined to be zero outside the sample and one inside, and its integral over all space is one). 
In the simulation, the CF is computed by direct pair counting (we have used the estimator 
introduced by fig]), and o~ 2 (r) by distributing N s spheres of varying radius with random 
centers and computing the number fluctuation. For the latter we have used N s = 2 • 10 4 
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spheres randomly distributed in the simulation volume and we define our estimator as 



where iVj(r) is the number of points contained in the i th sphere of radius r and (N(r)) is its 
average. One may see in Fig.| that for a CDM model, the cr 2 (r) is constant up to a scale 
r c « 0.06 (normalized to the simulation box L = 239.5Mpc/h), which is fixed by the turn-over 
scale of the PS H], and then it decays as r~ 4 . (The normalization has been performed by 
requiring that <j(r — 8Mpc/h) = 0.51 as in & we note however that a different normalization 
does not change our main results which concerns the functional behaviour of the real space 
properties with scale.) For the glass a 2 (r) ~ r~ A at any scale, as expected, while for the 
IC er 2 (r) decays as r -4 up to (A) and then it decays slower as r -16 : this change of slope is 
the effect of the displacements field. In Fig.|| we show the results for the two-point CF. For 
the CDM model used here || £(r) ~ const, up to r c and then after crossing zero, it goes as 
£( r ) ~ — r ~ 4 Hi- We find instead that £(r) for the IC and the glass are rather similar at 
small scales and they oscillate around zero: the peaks due to the first, second and third NN 
are clearly visible. Note that the average distance between NN is (A) = 0.003 (in normalized 
units) for both the glass-like distribution and the IC and it corresponds to the first peak of 
the conditional average density ~ (n)(l + £(r)) (see Insert Panel of Fig. 2). At large scales 
there is a slight difference between the two, which changes the behaviour of cr 2 (r). However 
in any range of scales there is no agreement between the £(r) and cx 2 (r) of CDM and IC. This 
is especially relevant for the nature of fluctuations in IC and CDM: for the latter one should 
see an over-density up to a scale r c followed by a peculiar under-density (with £(r) ~ — r~ 4 ) 
while for the first we find that they are a sequence of over-densities and under-densities which 
results in an oscillating £(r). 

Let us now discuss the second point of this paper. A lattice or a glass are gravitationally 
unstable with respect to small perturbations of their configuration. For this reason an in- 
teresting question concerns what drives the non-linear dynamics in these simulations. Once 
the Zeldovich displacements have been implemented, the new configuration does not have the 
perfect symmetry of the pre-initial distribution. Hence the question concerns: are discreteness 
effects dominant for the small-scales non-linear dynamics with respect to the force due to large 
scales smooth fluctuations ? In order to study this point, we have computed (see Fig.[|) the 
behaviour of the mean gravitational force on a particle due to all the particles contained in 
the sphere of radius r around, as a function of r. Firstly, we note that the force increases of 
a factor 10 approximately between the glass and the IC. In the case of the latter, one can see 
that the force due to the first shell is almost compensated by the second one and starts to grow 
at least until the sphere reaches the size of the box. Furthermore the force due to the first 
shells is one third of the force at large distance and fluctuations are of the same order than 
the average. One can therefore say that the contribution from the NN is not as important 
as in the Poisson case because the force doesn't reach its asymptotic value at scales ~ 2(A) 
[0- The average force increases with scale up to the box's size; however the fluctuations at 
all scales are large enough to conclude that small-scales discreteness effects are not negligible 
with respect to the large scale contribution. 

This last analysis would be totally useless without taking into account the initial velocity of 
the particles. Indeed if the particles were too fast they would not be trapped by nearby particles 
and discreteness would only make particles trajectories less regular. However, we expect the 
particles to be slow since we study CDM. In order to make this more quantitative, we did the 
following test: calculate the force on a particle due to particles in a sphere around, calculate 
the velocity of the particle in a direction perpendicular to the force and finally compare the 
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time needed to travel a certain distance with this velocity but without any force (t v ) and the 
time to travel the same distance with the force but without any initial velocity (tf). Our result 
is that the velocities are indeed small. For instance if we take (A)/ 10 for the distance to travel, 
one has that t v goes from ~ 3tf to ~ 10i/ when the radius of the sphere for the force goes from 
0.004 to 0.4. This shows that structures should be able to form even at the smallest scales 
r ~ (A) in the simulation. Indeed, at a distance (A)/10 from its initial position, the force 
from its NN could be already driving a particle, instead of its initial velocity. Discreteness 
effects are therefore important for what concerns the non-linear growth of structures at small 
scales and this is the relevant range of scales where one wants to study strongly or weakly 
(£( r ) ~ 0.1) non-linear regimes ||. (see |l6| for a more complete discussion of this point and 
of the whole time evolution of the simulations). 

Let us briefly summarize the discussion. The idea of the standard procedure to set-up IC 
in cosmological N-body simulations, is that one starts with a "uniform distribution" which, 
once discretised, can be a lattice or a glassy system. Then one gives correlated displacements 
to the particles in such a way that one gets a system which should behave like a continuous 
fluid with correlated density fluctuations which have the desired (i.e. representative of the 
continuous field) correlation properties. However, that one considers point distributions which 
introduce discreteness effects due to their intrinsic fluctuations and correlations. We have 
addressed the problem whether these discreteness effects due the imprint of the (correlated) 
fluctuations of the pre-initial point distribution are strong enough to super-seed the given 
correlations. Indeed, we have shown that this is the case and we have pointed out that 
the standard ad-hoc procedure used to create N-body IC does not give rise to the desired 
CDM-like correlation of density fluctuations. This implies that the small-scales non-linear 
dynamical evolution of the system is driven by fluctuations which arise from the particular 
ad-hoc procedure used to discretise the field. In this context, these fluctuations can be seen 
as finite size effects, and are completely different from CDM-like fluctuations. Moreover we 
have studied the behaviour of the gravitational force acting on an average particle with the 
result that the force due to the nearest particles is highly fluctuating and gives a contribution 
of the order of the large scale one. For this reason we conclude that discreteness effects can 
be significant relative to the smooth distribution and that they therefore play an important 
role in the small-scales non-linear evolution of cosmological N-body simulations. Finally we 
stress that most of procedures used to discretise a continuous field (e.g. threshold or random 
sampling) unavoidably introduces Poisson noise. In such a situation, due to the strong NN 
interactions, the large scale fluctuations play an even smaller role in the non-linear dynamics 
which would be dominated mainly by particle-particle interactions. 
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Fig. 1. - Behaviour of u 2 (r) for the pre- initial glass distribution, for the initial condition of the 
simulation (z — 50) and for the theoretical expectation for a CDM model. For the glass configuration 
<r 2 (r) ~ r , while for the IC it deviates from this behavior at r ~ (A) beyond which it behaves as 
r~ . For the CDM model instead cr 2 (r) ~ const, up to r ~ r c ~ 0.1 and then it decays as r~ 4 . (For 
the x-axis we have used normalized units to the box size L = 239.5Mpc/h). Note that there is no 
agreement at any scale, between IC and CDM. 



Fig. 2. - The behaviour of £(r) for the same distributions of the previous figure (for the x-axis we have 
used normalized units to the box size L — 239.5Mpc/h). For the glass and IC £(r) is oscillating around 
zero. For the theoretical expectation for a CDM one should find £(r) ~ const. > for r < r c ~ 0.1 
and then £(r) ~ — r~ 4 at larger scales. Note that the nature of the fluctuations in the IC and CDM 
distribution is rather different: for the first one has a sequence of over-densities and under-densities 
(i.e. an oscillating £(r)) while in the second case one would expect an over-density £(r) > followed 
by an under-density £(r) < 0. In the insert panel it is shown the behaviour of the conditional density 
~ [£(r) + 1] which makes clear the oscillatory nature of 



Fig. 3. - The behaviour of the average total force and its variance as a function of distance (for the 
x-axis we have used normalized units to the box size L = 239.5Mpc/h) for the pre- initial (glass) and 
initial particle distribution. For the glass configuration the force is very small and almost zero. For 
the IC, one may note that the NN gives a very fluctuating contribution to the total force of the order 
of the asymptotic one. Finally, note that the force does not converge at the box size. 
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